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ABSTRACT 

Investigating the evolution of disk galaxies and the dynamics of proto-stellar 
disks can involve the use of both a hydrodynamical and a Poisson solver. 
These systems are usually approximated as inhnitesimally thin disks using two- 
dimensional Cartesian or polar coordinates. In Cartesian coordinates, the calcu¬ 
lations of the hydrodynamics and self-gravitational forces are relatively straight¬ 
forward for attaining second order accuracy. However, in polar coordinates, a 
second order calculation of self-gravitational forces is required for matching the 
second order accuracy of hydrodynamical schemes. We present a direct algorithm 
for calculating self-gravitational forces with second order accuracy without artih- 
cial boundary conditions. The Poisson integral in polar coordinates is expressed 
in a convolution form and the corresponding numerical complexity is nearly lin¬ 
ear using a fast Fourier transform. Examples with analytic solutions are used to 
verify that the truncated error of this algorithm is of second order. The kernel 
integral around the singularity is applied to modify the particle method. The 
use of a softening length is avoided and the accuracy of the particle method is 
signihcantly improved. 

Subject headings: gravitation; methods: numerical 


1. Introduction 

Thin disks are common in the Universe as a result of the conservation of angular mo¬ 
mentum and efficient radiative cooling. The existence of central starburst rings (Lin et ah 
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2013; Seo & Kim 2014), bright and young stars formed along spiral arms (Elmegreen et al. 
2014), and substructures associated with bars and spirals (Kim et al. 2012; Lee & Shu 2012; 
Lee 2014) indicate that the self-gravity of gas is important to the evolution of disk galaxies. 
The formation of planets in the early phase of proto-stellar disks indicates that self-gravity 
of gaseous disks plays a role in shaping planetary systems (Zhang et al. 2008; Inutsuka et al. 
2010; Zhang et al. 2014). As a hrst approximation, these thin disks are usually studied using 
two-dimensional hydrodynamical simulations coupled with a Poisson solver. 

For a given mass distribution, the calculation of self-gravitational forces is a fundamen¬ 
tal, but a challenging aspect of computational astrophysics. For three-dimensional grid-based 
codes, several techniques have been proposed to improve the accuracy and the performance 
of calculations. Of the simplest ones is the use of the fast Fourier transform (FFT), which 
is suitable for both periodic and isolated boundary conditions (James 1977). The FFT 
techniques are fast, accurate and are suitable for both three- and two-dimensional calcula¬ 
tions. The multigrid relaxation methods are fast, flexible and have been used extensively 
when mesh rehnements are required (Hockney &: Eastwood 1988). However, the multigrid 
methods, which are by nature only for three-dimensional problems, cannot be reduced to 
two-dimensional calculations for an inhnitesimally thin disk as discussed in this paper. 

Compared to those techniques well developed for Cartesian coordinates, the calculation 
of the self-gravitational force in cylindrical and polar coordinates still requires further study. 
Yen et al. (2012) developed formulae for the calculation of these forces to 2nd-order accuracy 
for both Cartesian and polar coordinates. In this description, the Poisson integral is written 
in convolution form and is of linear complexity if the FFT is used. Unlike those Poisson 
integrals which are integrable in Cartesian coordinates, no closed forms were found for polar 
coordinates because the elliptic integral is involved. Consequently, while the 2nd-order ac¬ 
curacy can be achieved in Cartesian coordinates, the calculations in polar coordinates suffer 
from the presence of a singularity in the kernel integral, reducing the order of convergence 
to nearly hrst order. Convolution expressions of the Poisson integrals in polar coordinates 
is also adopted in Baruteau & Masset (2008, hereafter BM08) in their two dimensional hy¬ 
drodynamical calculations. However, in their formulae, the use of a softening length, though 
physically motivated, inhibits pursuing higher order accuracy. 

In this work, we develop a simple, but effective algorithm that increases the accuracy of 
the self-gravitational calculations in polar coordinates to 2nd-order. The proposed method 
retains a linear complexity since all the effort is directed to the preparation of accurate force 
kernels. The technique developed in this work is applied to improve the numerical accuracy 
of the particle method. The use of a softening length is avoided and the force kernel integrals 
in the neighborhood of a singularity signihcantly reduce the numerical error of the particle 
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method. 

This paper is organized as follows. The framework and assumptions adopted for this 
work are outlined in § 2. We develop the mathematical notations and formulae for the 
calculation of the self-gravitational force to 2nd-order accuracy and the modihed particle 
methods in both Cartesian and cylindrical coordinates, respectively in § 3. The 2nd-order 
method described for Cartesian coordinates is used in § 4, where two improvements for the 
evaluation of force kernels in cylindrical coordinates are elaborated. In § 5, detailed com¬ 
parisons between the numerical results and analytic solutions are discussed. We summarize 
our results and conclude in 5 6. 


2. Framework and assumptions 


The potential $ for a given distribution of gaseous density p in three-dimensional space 
satishes the Poisson equation below: 


V^$(a:;) = AttGp{x)^ 


( 1 ) 


where G is the gravitational constant and x denotes the position vector. Without loss 
of generality, we assume that G = 1 throughout this work. By imposing the boundary 
condition, 

lim $(a:;) = 0, (2) 

|a3|—>-oo 

the gravitational potential <h(a:;) can be cast in an integral form (Evans 1991; Binney & 
Tremaine 2008): 


^{x,y,z) 



IC{x — x,y — y, z — z)p{x, y, z)dxdydz, 


( 3 ) 


where {x,y,z) are Cartesian coordinates and JC = 1 /\/x^~^ry^^r^ is the kernel of the 
integral. In this paper, we restrict the discussion to the following expression for the density 
distribution: 

p{x) = a{x,y)6{z), (4) 

where 6 denotes the Dirac delta function and a{x, y) is the surface density dehned as: 


a{x,y) 


jp{x)dz. 


( 5 ) 


The integral form of the gravitational potential, i.e.. Equation (3), and the associated forces 
can be numerically evaluated through discretization in the computational domain. Yen et ah 
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(2012) have shown that uniform discretization in Cartesian coordinates and radially logarith¬ 
mic discretization in polar coordinates enable the self-gravitating forces to be expressed in a 
convolution form of a double summation. With the assumption that the density distribution 
is smooth, the linear approximation for the surface density in each cell increases the accu¬ 
racy of numerical solution. Using the convolution theorem (Bracewell 1999), a fast Fourier 
transform is applied to reduce the computational complexity from 0{N^) to 0(iV^ log 2 iV), 
where N is the number of zones in one direction. This method is a direct calculation of the 
self-gravitational forces, which is second order accurate in Cartesian coordinates, without 
necessarily invoking the use of artificial boundary conditions. 

One of the major advantages of the convolution integral approach is that most of the 
calculational effort is passed to the preparation of the kernels. That is, the more accurate 
the kernels, the more accurate the numerical solutions. We note that those kernels used to 
achieve higher order accuracy need to be prepared only once for a hxed grid and stored in the 
computer memory at the beginning of each simulation. In the following sections, we exploit 
the advantage of using the convolution integral and restrict the associated discussions to the 
plane of the disk, i.e., <F(a:,|/,0). Simple but effective approaches are proposed to improve 
the numerical accuracy and the order of convergence. 


3. A direct method of 2nd-order accuracy and a modified particle method 

In this section, we develop the mathematical notations that will be used throughout 
this work so that the material in this paper is self-contained. The expressions of formulae 
with 2nd-order accuracy are first derived in Cartesian and polar coordinates for readability 
and completeness. Based on the 2nd-order method, approximations are adopted to reduce 
the computational cost further by concentrating the mass of one cell at the cell center. The 
simplihed scheme is a modified particle method. The use of a softening length is avoided 
and the associated singularity problem is removed using the kernel integrals. Without in¬ 
creasing the computational cost, the modified particle-based method signihcantly improves 
the accuracy of the numerical solutions. 

In the following, we discuss in detail the calculations of the self-gravitating forces in 
the x-direction for Cartesian coordinates and in the r-direction for polar coordinates. In 
Appendix A, we provide the formulae for the calculations of the self-gravitational forces in 
the ^/-direction and 0-direction. The full expressions for the kernel integrals are also given. 
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3.1. Cartesian coordinates 


3.1.1. A direct method of 2nd-order accuracy 


Consider a calculational domain described hy D = [—M, M] x [—M, M] for some number 
M > 0, which is evenly subdivided with intervals in the x- and y-directions, respectively. 
Given a positive number iVd, we dehne Ax = 2M/Nd, Ay = 2M/Nd as the cell size in each 
direction and Xj+ 1/2 = —M + iAx, yj +112 = —M jAy as the cell boundaries, where i,j = 
0,iVd. The domain of each cell is then dehned to be Di j = [a:j_i/ 2 , Xi^i/ 2 ] x [|/j_i/ 2 , 2 /J+ 1 / 2 ] 
and the cell centers are Xi = {xi_i /2 + Xi+i/ 2 )/ 2 , yj = {yj_i /2 + yj+i/2)/2., with i,j = 1, ...,iVd. 
In total, the calculational domain is covered with iV| cells. 

The forces in the x-direction {Ffj) and ^/-directions {Ffj) are dehned at the center of 
cells and related to Equation (3) through the following relations: 


d 


Nd ATd 




^ i'=l j'=i 


i'=l f=i 
Nd Na 



d _ 

—lC{x-Xi,y-yj,Q)(T{x,y)<lx<ly ( 6 ) 

A /, ' 



D. 




lC{x - Xi,y - yj,Gi)a{x,y)<lx<ly. (7) 


The surface density a in cells, appearing in Equations (6) and (7) can be linearly approxi¬ 
mated by 

ct ( x , y) ^ + di, j,{x - Xi>) + 6f j,{y - yj>) ( 8 ) 

where cij/j', S^, ^ = da{xi',yji)/dx and S^, j, = da{xi',yji)/dy are constant in the cell Di'ji. 
With the linear approximation in surface density, Ffj with second order of accuracy can be 
approximated by (Yen et ah 2012): 


px ^ P^<^ p^,v 

hj ^ hj ^ hj ’ 

(9) 

Ni Vd 

i'=l j'=l 

(10) 

Vd Vd 

r^x,x _ \ ^ \ ^ cx ]^x,x 

i'=l j'=l 

(11) 

Nd A^d 

i^x,y \ ' \ ^ ry v^x,y 

i'=l j>=l 

(12) 


where 
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and 

(13) 

(14) 

(15) 

The first term in Equation (9) is the contribution if the mass enclosed within one cell were 
uniformly distributed and provides an accuracy of first order. On the other hand, the 
last two terms take into account the structure of the density distribution within the cell 
and, hence, provide an accuracy of second order. Equations (10) to (12) are convolution 
forms of double summations, which can be evaluated using EFT if the domain is uniformly 
discretized. Equations (13) to (15) can be integrated analytically and the detailed expressions 
are summarized in Appendix A. 


— i',j — j' 


i^x,x _ 









Ac,/ [(^ - + iy- 

{x - Xi){x - Xy) 

A,,/ [{x-Xi^ + iy- |/j)2]3/2 
{x - Xi){y - yf) 

A/ [(^ “ + iy- 


dxdy, 

dxdy, 

dxdy. 


3.1.2. Modified particle method 


The 2nd-order scheme described above involves double summations with three types of 
force kernels, i.e., and when calculating forces in x-direction. 

The computational cost can be considerably reduced if a further approximation is adopted 
for Equation (6): 


F: 




(Vd 

+ (%' - J Jd^,., 

Ni TVd 


Xjf Xi 



a{x,y)dxdy, 


i'=l j'=l 


(16) 

(17) 


where e > 0 denotes the softening length. Equation (17) is equivalent to placing a particle 
with mass Mj/j/ = f cr(x, y)dxdy at the cell center (x^/, yp) and the Ff ■ is approximated 
as a result of a direct summation. That is, the A^-body calculation uses Mj/j/ = ai>j>AxAy. 
We note that the expression of Equation (17) is still a convolution form of double summation 
with a kernel = {xy —Xi)/[e^ + (x,/ —Xifi + {yj' — yjY]^^'^- A EFT can be applied to 

reduce the computational cost compared to the use of direct summation. A nonzero softening 
length e is usually introduced in the denominator of j_j, to avoid the singularity when 
x' = X. This calculation involves only one double summation, but at the expense of an order 
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of accuracy. One may expect the calculation of Equation (17) is roughly three times faster 
as compared to that of Equation (9). 


Introducing a softening length is not a desirable feature since forces are distorted, which 
reduce the accuracy of numerical solutions. However, the problem of a singularity does 
not exist in Equation (13) since it is integrable. This suggests a way to avoid the use of 
a softening length and an improvement to the kernel function is possible. The 

improved force calculation T)®- is proposed as follows: 


^ a- ^ M-, ■, 


Nd Nd 

EE 

i'=l j'=l 


iVd A'd 

i'=l j'=l 


(18) 

(19) 


where 


j^X,p 


T^x,c 


[{Xi' 

0 , 


Xi' Xi 

XiY + (%■' - 2/i)^]^/^’ 


rr. .r^’O I ca: ^X,x , ry J^x,y 


OT j ^ f 

otherwise. 


( 20 ) 

( 21 ) 


px,corr ^ correction term that takes into account the gravitational force contributed from 
Dij, which is set to be zero in Equation (20), reflecting the idea that a particle does not 
feel its own gravity. On the right hand side of Equation (21), the hrst term is zero due to 
the symmetry of the cell, the last term is also zero since the integrand in Equation (15) is 
an odd function of y with respect to the cell center. Only the second term that involves the 
gradient of the surface density in the x-direction is in general nonzero. Compared to the 
double summation, the cost of calculating jg very small, but it signihcantly improves 

the accuracy of the numerical solution as will be shown in Section 5. The particle method 
that does not suffer from the problem of singularity and involves higher order correction for 
x' = X is called the modihed particle method in this paper. 


3.2. Polar coordinates 


3.2.1. A direct method of 2nd-order accuracy 


Corresponding to Equation (3) for Cartesian coordinates, the potential function in polar 
coordinates can be expressed as: 


$(r,0,2) = - 



/C(r, r,(j) — (f), z — z)p{r, z)rdrd(j)dz, 


( 22 ) 
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where (r, 0, z) are polar coordinates and /C(r, r, </>—0, 2 ;— 2 ;) = 1/ ^/r'^ — 2rr cos(0 — 0) + + (^r — 2 ;)^. 


Consider the compntational domain described hy 71 = 7Z^ U 7?.® U 71 j, which is the 
nnion of three parts 7Z‘^, 7Z^ and TZj. Here, TZ'^ represents the destination domain, where 
the resnlting gravitational forces inclnde the contribntions from the whole compntational 
domain 7Z. 7Z^ is the sonrce domain, in which the mass gravitationally inflnences R^. TZj 
contribntes the gravitational forces associated with the origin of the calcnlational domain 
and its surronndings, which is not inclnded in 7Z'^ and 7Z^. The discretization of the domains 
7^'^, 7Z^ and TZj adopted in this work is described below. 


7Z^ = [Md , X [0, 27r] for some nnmber > Md > 0. The radial direction is 
discretized in logarithmic form and the azimnthal direction is evenly subdivided. Namely, for 
a positive integer we dehne A0 = 27r/N'd, (3 = (M^m/Md)^/^'!, rj+ 1/2 = 0j+i/2 = 

jA(j), i,j = 0,...,iVd, Vi = (ri_i /2 + rj+i/ 2)/2 and (j)j = {^j- 1/2 + 0j+i/2)/2, where = 
1, ..., Ad. The destination domain is covered with N\ cells dehned by TZf j = [nj_i/ 2 ,'rj_|_i/ 2 ] x 
[0j_i/2, 0J+1/2] for i,j = l,...,Ad. We note that the arrangement of TZfj does not cover 
the region r < Md. Extra cells that cover r < Md should be included in the region 
7Z^ = [Mj^, Md] X [0, 27r]. The region 7Z^ is discretized in the same way used for discretizing 
TZ’^, i.e., using the same [3 and A0. Without loss of generality, the region 7Z^ is discretized 
with cells, and with Mf^ = (3~^’‘M-^. We dehne r'j_,_i /2 = /7*Md, for i = —As, ...,0, and 
Tj = (rj_i /2 +?"i+i/2)/2 for i = — Ag + 1,...,0. Since the discretization in the azimuthal 
direction is directly inherited from that used for TZ'^, no special care is required. The cells 
dehned by TZlj = [rj_i/ 2 , ri+ 1 / 2 ] x [0j_i/2, 0j+i/2] for i = -Ag + 1,...,0 and j = l,...,Ad 
are used to cover 7Z^. Finally, cells, TZj = [0, Mf^^j x [0j_i/2, 0j+i/2], should be included to 
take into account the contribution from the vicinity around the origin. For simplihcation of 
notation, we denote TZ^p^^j = TZj and TZij = [ri-i/ 2 Ri+i/‘ 2 \ x [0j_i/2, 0j+i/2] for the ranges 
of indices i = —Ag + 1,..., Ad and j = 1,..., Ad. 

The forces in the r-direction {F^j) and 0-direction {Pfj) are dehned at the center of 
cells (r,, 0j) G TZ'^ and related to Equation (22) through the following relations: 


— 




^/C(r, Ti, 0 - 0j, 0)(T(r, 0)rdrd0, (23) 


7C(r, Tj, 0 — (pj, 0)cr(r, 0)rdrd0, (24) 


where the surface density a{r, 0) in TZi^j/ is linearly approximated by: 



a(r,0) 




(25) 
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where SI, j, = da{ri/, (j)j')/dr and Sf, = da{rii,(j)j/)/d(j) are constant in the cell TZi'ji. 
With the linear approximation in the surface density, F[j with accuracy of 2nd-order can be 
approximated by (Yen et ah 2012): 


F 




j!r,0 


^ 1,7 ^ 1,7 ; 


(26) 


where 








and 


Nd 

Nd 


0 Nd 

Nd 



+ 

1 


' + CT-ATs 


i'= 

i'=i 

1 

'=-7Vs+li'=l 

i'=i 



■A'd Nd 


0 Nd 

Nd 


F 

EE 

cr ^r,r 

1 \ ^ \ ^ rr ir’7',7’ 

T / ^ 

.-k + E^- 

Ns.r^i+N,,j- 


j'=i j'=i 


i' = — Na + l j' = l 

i'=i 


Nd 

Nd 


0 A^d 

Nd 


EE 7, 

+ 

1 

E 

+ E ^-Ns,3 

''^i+NaJ-j'^ 

i^= 

j'=l 

i 

= -Na + l j' = l 

j'=l 



(27) 

(28) 
(29) 


v^r,0 


i^r,r 


^r,4> 



f(r,-fcos(^-</>,)) - 

[^2 F rj — 2fri cos{(j) — ’ 

f(rj-fcos(0-0j))(f-rj/) - 

rj[f2 + r'f — 2fri cos(0 — 

fjn - f cos(0 - - (pf) 

[^2 _|_ ^2 _ cos (0 — 


(30) 

(31) 

(32) 


On the right hand side of Equations (27) to (29), the hrst terms are the self-gravitating 
terms, the second terms are the “gravitational interaction terms” contributed from 77® and 
the last terms are the contributions from the vicinity around the origin TZj. Although the 
hrst two terms can be mathematically combined into a single double summation, in practice, 
we treat those three terms separately. The self-gravitating terms involve integration around 
the singularity, while the interaction terms do not. The third terms, though in a convolution 
form of a single summation, are not compatible with the hrst two terms and therefore treated 
separately as well. EFT can be applied to all these calculations to reduce the computational 
cost and keep the complexity at 0{N‘^ log 2 N). 


Some properties regarding Equations (27) to (29) are worth mentioning. First, the order 
of accuracy relies on how accurate the force kernels, i.e.. Equations (30) to (32), are evaluated. 
As pointed out in Yen et al. (2012), no closed forms are found for the and 

j-j" The integration involves an elliptic integral that can only be evaluated numerically. 
Moreover, the presence of a singularity function in terms of ln(l — cos(0)) degrades the order 
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of accuracy to first order. It is therefore desirable to improve the accuracy of the force kernels 
given the fact that Equation (26) involves three integrals, which are originally dedicated to 
reach an accuracy of 2nd-order. These considerations will be addressed in the next section. 
Second, Equations (30) to (32) can be used to evaluate “gravitational interaction” between 
two grid patches, i.e., the destination patch and the source patch 7?.®, which are partially 
or completely separated. As long as the their domain discretization shares the same f3 and 
A0, EFT can be applied to reduce the computational cost. In the case that and 7Z^ 
are completely separated, we do not need to worry about the singularity associated with 
ln(l — cos(0)). We note that a constant spatial shift between two patches in r and (j) is 
allowed, since the spatial shift only contributes constant phase shifts in the Fourier domain. 
Third, in the case that TZij = TZi/ji, where the singularity occurs, the values of /Cpo, ICq’q 
and /Cg’o are invariant. For example, casting Equation (31) in the form 


/C 


r,r _ 

0,0 ~ 


rA(t,/2 rW(/3+l) 

J-A ^/2 A/C/^+i) [v'^ + l- 2r/ cos(0]^/^ 


(33) 


where rj = f/vi and ^ = (j) — (pj, /Cq’q is a constant for all i' = i and f = j since A0 and 
P are constants. This is a useful property given that the cell sizes are not uniform in polar 
coordinates. This indicates that one can place the effort on evaluating the elliptical integral 
for one specihc cell which contains a singularity and apply the result to all other cells that 
also contain a singularity. 


3.2.2. Modified Particle Method 


The modihed particle-based method can also be applied to polar coordinates. Similar to 
Equation (16), we can approximate Equation (23) to further reduce the computational cost 
at the expense of the order of convergence. Corresponding to Equation (18) for Cartesian 
coordinates, the approximation in polar coordinates is written as: 


Fii « 'y,.iril + TiSliK.Z, + stP«t+ E 


Na Na 

E E 

i'=-7Vsi'=l 


Na Na 

i'=-Ns j'=l 


(34) 

(35) 
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where 




E 





Vji cos{(j)j> - 0j)) 
2riiri cos{(j)f - 


+ ns: 


^r,r 
hJ^O,0 


X<f> irr,<l> 

hj^0,0 • 


i^i' OT j ^ f 

otherwise. 


(36) 

(37) 


Here, denotes the gravitational force contributed from T^ij, which is set to zero in 

Equation (36). It can be also shown that /Cq’q = 0 since the integrand in Equation (32) is 
an odd function of 0 with respect to the cell center. The hrst two terms on the right hand 
side of Equation (37) are in general nonzero. Since the shape of cells in polar coordinates 
is not symmetric with respect to the center of cells, involving the correction term cTj^j/CQ q is 
particularly relevant in polar coordinates. The cost of computing is small compared 

to the calculation of a double summation. 


4. Singular integration method (SIM) 


The mathematical formulas developed for 2nd-order convergence have been shown in 
Section 3.2. The lack of closed forms for Equations (30) to (32) dictates that the order of 
convergence relies on the numerical methods used for the integrations. Yen et ah (2012) 
evaluate the force kernels using the trapezoidal rule with one trapezoid as follows 








-HI 

-(4> - 4>f)H\ 




V-1/2 ’ 


^i' + l/2 i‘?^j' + l/2 _ Ifrfi 

V - 1/2 1/2 ~ 



^i' + l/2 Yj/+l/2 

V - 1/2 4 j /_ i /2 ’ 


(38) 

(39) 

(40) 


where the notation /(■)]^ = (/(&) + /(a))(fe — a)/2, and the exact expressions of 77^ and 
are given in Appendix A. It is natural to utilize more than one trapezoid to improve the 
accuracy. Specihcally, 


-^tpz 





+1/2 

-1/2 

Em+\ 

4k ’ 

(41) 


m=l 






A^tpz _ 



A' j^r,0 

4k 


^r,r 



+1/2 

-1/2 

(42) 


m=l 



* 



^tpz _ 







0 

— (pj 

V+l/2 l^k+l 

J V-1/2 Jflik ’ 

(43) 


m=l 
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where = 0j'-i/2 + (m - 1)A9, A9 = {(j)j'+i /2 - 0j'-i/2)/^tpz and Ntp^ > 1 denotes the 
number of trapezoid used for the evaluation. As will be shown in Section 5, this consideration 
signihcantly improves the accuracy of numerical solutions in the self-gravitating case. We 
emphasize that evaluating Equations (41) to (43) using Wpz > 1 does not increase the 
computational complexity of the method, since those kernels are calculated only once at the 
beginning of simulations. 

Another issue associated with Equations (30) to (32) is that when i = i' and j = j', a 
singularity in the form of ln(l — cos(0)) is involved in those integrals. In this situation, forces 
evaluated using Equations (38) to (40) can have incorrect signs, which not only degrades the 
order of convergence to hrst order but also deteriorates the accuracy of numerical solutions. 
Thus, special care is required for the evaluation of JCqq, and /Cqq. Fortunately, Equa¬ 
tion (33) suggests that special care need only to be taken once for one specihc cell and the 
result can be applied to other cells. 

As shown in Figure 1(a), we cover a specihc fan-shaped cell using Cartesian cells. The 
fan-shaped cell can be characterized by: 


Ax = 
Ay = 



cos 



(44) 

(45) 


where (r^, 0m = 0) denotes the cell center and (Ar, A0) dehnes the size of the fan-shaped 
cell. Since (Ar, A0) decreases with increasing Ad, the number of Cartesian cells used to 
cover the fan-shaped area should increase accordingly to well resolve the curved fan-shape. 
A good rule of thumb is to cover Ax and Ay with roughly 10 Cartesian cells in x and y 
direction, respectively. As a result, the corresponding number of Cartesian cells used to cover 
the fan-shaped cell in this work is lOAd in the x direction and 5Ad in the y direction. The 
surface density of those Cartesian cells that lie outside the fan-shaped area is set to zero. 

We evaluate the forces at the cell center, (rm,0m), associated with (a) uniform surface 
density (j(x, y) = 1, (b) a{x, y) = f—r^ and (c) (j(x, y) = 0—0m- Since we are only interested 
in the self-gravitating forces at a specihc point, (rm, 0m), evaluation using Equation (9) that 
involves only three double summations would suffice. The force calculated for Figure 1(a) 
corresponds to /Cqq due to unit surface density. Figure 1(b) corresponds to rmJCl’o due to 
the unit radial slope in surface density, and Figure 1 (c) corresponds to /Cqq = 0 due to the 
unit azimuthal slope in surface density. Figure 1(c) is more relevant to /Cq’q as shown in 
Appendix A. We note that when applying the force, associated with Figure 1(b) to 
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other cells, a factor ri/vm is required, i.e.. 


r,r ^ _;_2_ pr,r 

'^tj ry^ T'mi4>n 

• m 


p^r.r 


(46) 


due to the factor appearing in front of the square bracket in Equation (28). Hereafter, we 
call the algorithm described in this section as the singularity integration method or SIM in 
short. 


5. Results 


In this section, we verify the accuracy and the order of convergence proposed in this 
work by comparing the numerical solutions with examples that have an analytic solutions. 
For Cartesian coordinates, we focus on the accuracy improvement for the modihed particle 
method, while we show improvements in both the accuracy and the order of convergence for 
polar coordinates. 

We investigate the numerical error in the destination domain TZ'^ and V using the 
following dehnitions of error: 


L°° 


Nd Nd 


jpEEk. 

d j=\ 


num _ pani 


EEkr--f;Tr ■ 

i=i j=i J 

|) for i, j G or D, 


1 

M- 

d 1=1 j= 


max( 


(47) 

(48) 

(49) 


where L}, are the one norm, two norm and maximum norm of error, and E-™™, 

are numerical and exact forces at locations indexed by (f, j), respectively. When using L} 
and L^, we evaluate the total variation and energy in a global sense, while using we focus 
on the convergence of maximum error in a pointwise sense. 


5.1. Examples with analytic solutions 

Direct comparisons between numerical with analytical solutions are desirable to demon¬ 
strate the effectiveness of a numerical method. In this work, we are concerned with the 
accuracy and the order of convergence of the proposed algorithms. For these purposes, the 
selected disk models need to fulhll the following criteria. First, the disk is supposed to be 


14 


infinitesimally thin and has closed-form solntions for the self-gravitational forces and the 
potential. Second, the size of the disk should be hnite to be fully enclosed by the hnite 
calculation domain. Third, the mathematical form of the density distribution should be 
sufficiently smooth, i.e., higher order derivatives behave well in the calculation domain, for 
analyzing the order of accuracy. Disk models that do not fulhll these criteria are not useful 
in understanding the properties of the proposed algorithms. 

Only few inhnitesimally thin disks are found to have corresponding closed-form solutions 
of self-gravitational forces, e.g., the Mestel disks (Mestel 1963), the exponential disks and the 
generalized Maclaurin disks (Schulz 2009). Among them, only the density-potential pairs 
discussed by Schulz (2009) satisfy the first two criteria above mentioned. Schulz (2009) found 
closed-form solutions in cylindrical coordinates for the hrst three members n = 0,1, 2 of the 
family of hnite disks with surface density, described by: 





for r < a 
for r > a, 


(50) 


where r = + y'^ and a is a given constant describing the size of disks. It can be shown 

that even the disk, which is smoothest among the three, has a singularity in the second 
derivative of Equation (50) at r = a. In other words, disk is not sufficiently smooth 
along the disk edge and will degrade the order of convergence in terms of the maximum 
norm error. To circumvent this issue, we generalize the closed-form solutions for arbitrary 
positive integer n > 0 as the following. 


In general, for a given integer n > 0, the surface density of ctd„+i can be associated with 
though the following recursive relation; 


az)„+i(r;a) 


2n -|- 1 

Q,2n+1 



«)dd. 


(51) 


where a serves as a dummy variable for the integration. Due to the linearity of the Poisson 
equation, the corresponding radial force at the mid-plane has a similar recursive relation: 




2n + 1 

Q,2n+1 



^ In 77 r,ana/ ^ \ a ^ 

“ (’■io)do. 


(52) 


Without loss of generality for the discussion in this work, we assume (Jq = 1 and G = 1. 
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The closed-form of the radial force has the following general form: 


/r N 2n-l 

Y Va 


n—1 


a) = 


/ f \2n—l 

n—1 


k=0 

k=0 

E (7) Vl - Hry- 


for r < a 


\ k=0 


for r > a 


(53) 


where T 2 k is the Chebyshev polynomial of the hrst kind of order 2k and (a^fc+n ^ 2 k) ^^e 

coefficients associated with the odd and the even order of Chebyshev polynomials, respec¬ 
tively. The coefficients have the following recursive relation with (a"’, 6 "’): 


and 


a^++\ = ( 2 n + l)^, 


..n+l 


tn 


= ( 2 n + l) + ^ , 


3 = 

J = 0 




a. 


'2k+l 


lE 


a, 


n 


2k-l 


4:^ 2k+ 2 4^ 2k 

k=j k=j 

n—3 n -t n—1 


1 ^ ®2fc+5 


4^ 2k 

k=j 


E 


^2k 


2 ' 2 ^ (2A: + 2)2 

k=j 


^ n—2 2 

_£ ^2k+3 

2 4^ 2k+ 2 

k=j 

-I n—3 in 
_ £ ^2fc+4 

2^^(2k+ir 


(54) 


j = 0 , ...,n 


'' 6 ;+'= ( 2 n + l) 5 , j = 2 ,...,n 

b"*' = ( 2 n + 1 ) £ + |£ j = 1 

6r‘ = (2n + l)(^-££ j=0 

3 = 


1 / 5” 
^ = 1 / 


tn 

^2j+2 


(55) 


J 


with a{ = 1 and 69 = —1. The derivation of Equations (54), (55) are detailed in Appendix 
B. By using a ao^. disk, the surface density has smoothed n — 1 order of derivative at 
the edge of the disk. In the following, we adopt n = 2 and 5 in illustrating the issue 
associated with the smoothness of the surface density, i.e., the third criterion, and justify 
that the SIM is of nearly 2nd-order accuracy. Since the Poisson equation is linear and 
those examples considered in this paper involves all Fourier modes in both the radial and 
azimuthal directions, the conclusions drawn from this work are general and applicable to 
any other smooth density distribution. 
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5.2. Results of Cartesian coordinates 

The comparisons of the calculated radial force in Cartesian coordinates for different 
methods are shown in Figure 2 for udj with a = 0.25 using N = 128. In Figure 2(a), the 2nd- 
order method and the modified particle method (denoted as particle-|-) have better numerical 
accuracy compared to that of particle method without correction (denoted as particle). The 
absolute value of the radial force is significantly underestimated in the last case, which 
neglects the density gradient in one cell. Since the gradient is negative, i.e., density is higher 
toward the center of the disk, one may expect a radially inward force contributed from 
inside a cell. The improvement is best shown in Figure 2(b) using relative error defined by 
|^r,num _ ^r,ana | ^ | ^r’,ana| _ Compared to the particle method without correction, the absolute 
error is reduced by an order of magnitude if the density slope within the cell is taken into 
account, and an additional factor of five improvement is found in the 2nd-order method. 
This is significant due to the inverse square law of the gravitational force. 

Figure 3 shows the one norm error of radial forces {Ll) as a function of cell number N^. 
The decrement in L^. with increasing indicates the convergence of all methods. The slope 
corresponds to the order of convergence as indicated by the solid line for the Ist-order and 
the dashed line for the 2nd-order. This figure shows that, in general, the 2nd-order scheme 
is the most accurate method compared to others and indeed has numerical convergence of 
nearly 2nd order. Although both the particle-based methods have numerical convergence of 
nearly first order, the numerical accuracy of modified particle method is better than that of 
the particle method by one order of magnitude. 


5.3. Results of polar coordinates 

We have three cases for cylindrical coordinates. In the first case, we demonstrate the 
order of convergence of the methods in the absence of a singularity. To do so, a disk with 
a = 0.006 centered at the origin is employed in the source domain 10“^] x [0,27r] 

and TZ = [0, Mf^] x [0, 27r] that encloses the origin, keeping the destination domain TZ’^ = 
[10“^, 1] X [0, 27r] devoid of mass. We adopt = lO”"^ when = 8 and the value of 
consecutively shrinks to roughly half the size whenever Ah is doubled. 

Figure 4 shows the one norm and the maximum norm errors of the radial forces as 
a function of N^. The convergence of both the particle method and the original method 
proposed by Yen et ah (2012) are of nearly 2nd-order. The 2nd-order convergence of the 
particle method is general if the destination domain TZ'^ is devoid of mass to avoid the 
singularity involved in the integrals Equations (30) to (32). This can be understood as 


17 


described below. 

Consider the location {ri,(j)j) G TZ'^ where one feels the gravitational force in the radial 
direction from 77®, 

* iJ 


pr 



/ , r'[rj — r'cos(0'— 0i)l , ,,,, 




[r'2 + rf — 2r'ri cos(0' — 
substituting the following relations into Equation (56): 

r' = Tj/ + r, 

0 ' = (py + 0 ^ 

cr(r', 0') ^ CTi/j/ + 5[/j/(r' - r^/) + Sf,j,{4) - pf) 

= aiij> + 6l,j,f + 5)^ 

cos(0' — (pj) = cos(0 + (pj! — (pj) = COs{(pf — (pj) 

^3 


+ 0 ( 0 ^ 


sin( 0 j/ — (pj 


6 


+ 0 ( 0 ® 


(56) 


(57) 

(58) 

(59) 


(60) 


It can be shown that: 




-ri'[ri - Ti! cos(0j/ - (pj)] 


A0y/2 Ar.,/2 


[rl + r‘f - 2riiri cos{(pj' - 0j)]^O J J 

-A0y/2 -Ari,/2 

+ 720 + 73!"^ + 740^ + 75^0 + O(f0^))dfd0 

= + 0 ((Ar.)^ + (A,#,,,)^)) 

[rji + rt — 2rj,rj cos(0j/— 0j)J+"^ 


(cTi/j/ + + 5)^_^.,0)(1 + 7 ir 

(61) 


+5[,,^.,Ari,A0,v(O((Ar,O2 + (A0,O2))) + sf,^^,Ar,AP)j,{0{{Ar,)^ + (A0,v)2))) 
-[ri - Ti! cos{(pj> - (pj)] 


[rl + rf - 2ri>ri cos{(pj' - (pj)]^/"^ 


Miiji + 0 ((Arj/) + (A 0 j/) )), 


(62) 


where Mi'j' = aiijirifArifA(pji. Equation (62) indicates that the accuracy of the particle- 
based method is in general of 2nd-order in the absence of a singularity. Equation (61) involves 
the Taylor expansion of the denominator in Equation (56). When approaching Equation (61) 
from Equation (56), we have to assume that the distance \x — x'\ ^ Ar*/ in order to have a 
reasonable speed of convergence. This assumption breaks down in the self-gravitating case 
that involves an integration around a singularity that reduces the order of accuracy. Figure 4 
also shows that the particle method seems to be more accurate than the original method 
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proposed in Yen et al. (2012). The difference comes from the use of central difference at the 
disk edge. The disk mass does not vanish to zero at r = a in the latter case. In spite of this, 
the error is reduced commensurate with that of 2nd-order. 

In Figures 5 and 6 we show the maximum norm errors for the ctdj disk and the (Td^ 
disk, which are centered at [xc = 0.5, |/c = 0.1) and with a = 0.25. In contrast from 
the first case, the center of the disk is shifted from the origin of cylindrical coordinates, 
providing nonzero azimuthal forces in the domain TZ'^. Both cases involves the integration 
around singularities and is therefore useful for testing the algorithm discussed in Section 4. 
Since L°° is particularly helpful for monitoring the convergence of the maximum error in the 
computational domain, only L°° is shown in the form of figure. The order of convergence of 
different algorithms associated with L^, and L°° are tabulated in Table 1 for both 
and CT/jg disks. 

The top panel of Figures 5 and 6 shows the maximum norm errors of the radial forces, 
while the bottom panel shows that of the azimuthal forces. Four different numerical algo¬ 
rithms are shown in the plots. The open circles are the results obtained from the particle 
method, diamonds are from the SIM described in Section 4 using iV^pz = 19, the asterisks 
are from the modified particle method and the triangles are obtained from the algorithms 
described in BM08. When implementing the last algorithm, a softening length e = 0.015r 
is adopted as that was used in BM08 for the model with = 64. For a fair comparison, 
the size of softening length used for other models is scaled linearly with the mesh size. For 
instance, the softening length used for Nd = 128 is e = 7.5 x 10“^r, while for Nd = 32 is 
e = 0.03r. We note that the use of a softening length in BM08 is well guided by the physical 
consideration of disk thickness (Muller et al. 2012). These plots and Table 1 show that both 
the particle method (denoted as particle) and the method used in BM08 have genuine first 
order of accuracy, since the error decreases linearly with the cell size. Compared to the re¬ 
sults of the particle method and BM08, the modified particle method (denoted as particle-|-) 
effectively reduces the numerical errors. The difference between the particle, the BM08 and 
the modified particle methods is only on the treatment of a singularity involved in Equa¬ 
tions (30) to (32). Although the improvement on the numerical accuracy is significant, the 
modified particle method is still of first order convergence. Another significant improvement 
is achieved when the SIM described in Section 4 is adopted. In addition to the treatment 
of the singularity, numerical integration of Equations (30)-(32) is improved using the trape¬ 
zoidal rule with more than one trapezoid. In this work, we find that using Wpz > 5 gives 
reasonable numerical accuracy and the results will not signihcantly change with increasing 
Wpz. In the case with a/jg, the order of convergence is nearly second order. For the model 

with Nd = 1024, SIM is two order of magnitude more accurate than the particle and 
the BM08 methods, and one order of magnitude more accurate than the modified particle 
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method. Meanwhile, the computational complexities of all these four methods are the same 
0{NHog,N). 

From Figure 5 and Table 1, for the disk, SIM seems to have only roughly 1.5 order 
of convergence in terms of L°° while it has nearly 2nd-order accuracy in terms of global error 
measurement using and L^. This indicates that the maximum error in a £)2 disk converges 
slower than that in In Figure 7, the maps of absolute error are used to investigate the 
distribution of error for and disks. These maps are produced using the SIM. The 
values of error are color coded as indicated by the corresponding color bars. It is evident that 
the major errors are concentrated at the edge of the aD 2 disk. On the other hand, the errors of 
the (jDg disk are smoothly distributed over the disk. The reason is that the aD 2 disk described 
by Equation (50) cannot be well approximated by the linear expansion Equation (25) when 
approaching the edge of the disk. That is, for n = 2, a singularity develops in the second 
derivative of Equation (50) at r = a. This phenomenon does not occur in the (Td^ disk since 
it has a smooth second derivative throughout the computational domain. Figure 7 justihes 
the third requirement for the disk model as mentioned at the beginning of Section 5.1, i.e., 
the density distribution of the model disk needs to be sufficiently smooth for analyzing the 
order of accuracy. 


6. Discussion and summary 

Equation (62) indicates that the particle method is of nearly 2nd-order convergence 
in cylindrical coordinates when \x — x'\ 3> Ar*/. Similar argument and conclusion can 
be also applied to Cartesian coordinates. This is a desirable property for the calculation 
of gravitational interaction between patches if their domains are mutually exclusive and 
separated. This situation is commonly seen in a numerical code featured with adaptive 
mesh rehnement. We suggest to use the particle method for the gravitational interaction 
between two separated patches, and apply the SIM only for the self-gravitational forces inside 
a patch. In the following, we summarize this work. 

Building on the work of Yen et al. (2012), the self-gravitational force calculation for an 
inhnitesimally thin gaseous disk in cylindrical coordinates is improved. The original method 
proposed for the cylindrical coordinates (Yen et al. 2012) is only a scheme of approximate 
hrst order in convergence. We identify two sources of error that degrade the numerical ac¬ 
curacy and the order of convergence. One arises from the use of the trapezoidal rule for the 
kernel integration with only one trapezoid, the other is due to the presence of singularity in 
the kernel integral when x' = x. The former issue is resolved by increasing the number of 
trapezoids for the integration, while for the latter we adopt the 2nd-order method in Carte- 
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sian coordinates, which is free of singularity, to evaluate the integral around a singularity in 
polar coordinates. We prove that the result of integration obtained for a specific cell can be 
applied to other cells if the radial direction is discretized logarithmically and the azimuthal 
direction is discretized evenly. These two improvements significantly reduce the numerical 
error and result in nearly 2nd order convergence. 

A similar consideration is applied to the particle method. We show that the particle 
method is of 2nd-order convergence in the absence of a singularity. When singularities are 
involved in the self-gravitational calculation, the use of softening length degrades the nu¬ 
merical accuracy. Thus, we propose to incorporate the force integration around a singularity 
as for the SIM to the particle method. As a result, the accuracy is significantly improved 
for (J 02 and disks. However, this correction is not sufficient for improving the order 
of convergence since neglecting the detailed distribution of mass in the surroundings of a 
singularity introduces an error of hrst order. 

These considerations do not increase the computational complexity 0(iV^ log 2 A^) since 
all the efforts are focused on improving the accuracy of force kernels, which are only cal¬ 
culated once at the beginning of the simulations. The method for self-gravitational forces 
presented here can be similarly applied to the gravitational potential (Yen 2014). 

Figures 5 and 6 show a different rate of convergence in terms of L°° for the ao 2 and 
aDs disk, respectively. By noting that the maximum errors are concentrated at the edge of 
the disk as shown in Figure 7, we conclude that the rate of convergence is related to the 
smoothness of the mass distribution. We note that the second derivative of aD 2 disk does not 
exist at the edge of disk. This can be also understood mathematically from Equations (23) 
and (24). The integrations involve two parts, one is associated with the force kernels and the 
other is associated with the surface density. Numerical experiment shows that the results 
for ct £)2 disk do not significantly change as the number of increases from 19 to 39. This 
indicates that an improvement on the force kernel integral by increasing the number of 
cannot increase the numerical accuracy any further. Thus, the lack of 2nd order behavior 
in terms of L°° as shown in Figure 5 can only be associated with the term related to the 
approximation of the surface density, i.e.. Equation (25). When applying Equation (25), 
we implicitly assume the underlying density is sufficiently smooth so that the error of the 
approximation is of 2nd order. However, this statement is not true at the edge of the aD 2 
disk. Thus, we should not expect a 2nd-order accuracy in terms of for disk. The 
smoothness assumption may seem to be a limitation of the SIM. However, given the fact that 
those density values are given for a set of discretized points, without a priori knowledge about 
the functional form, the linear approximation using Equation (25) is perhaps a good way to 
reach second order accuracy. We note that the smooth regions can be approximated with 
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higher order accuracy, while those regions with discontinuities can be only approximated 
with lower order accuracy to avoid the Gibbs phenomenon. For instance, a slope limiter is 
designed to reach second order accuracy in smooth regions and to avoid numerical oscillations 
around discontinuities when solving hydrodynamic equations with the Godunov method. 

We have shown that the use of a softening length reduces the accuracy of a Poisson solver 
to first order (Baruteau & Masset 2008; Li et ah 2009). To be commensurate with the second 
or higher order accuracy of hydrodynamical solvers, the use of a softening length should be 
avoided. One limitation of the method described in this work is the use of logarithmic radial 
grid. This grid configuration may be useful if the computational domain requires a large 
spatial range, e.g., a protoplanetary disk. Readers who are interested in using a uniform 
discretization in the radial grid should refer to the work by Li et ah (2009). 
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A. Full expressions of kernels 
A.l. Cartesian coordinates 


The calculation for Ff- is fully analogous to that for With the linear approximation 
in surface density, F^ - can be approximated by: 


rv ~ pvfi T?y^y i 


(Al) 
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where 


Ni A^d 

i'=l j'=l 
Nd A^d 

T^y,y _ xy )ry,y 

^i,j 2-^ 2-^^i',j'^i-i',3-3"’ 
i'=l j'=l 
Nd ATd 

r^y,x _ XX ^y,x 

^i,3 - / . / ."i',3''^i-i',i-i'^ 
i'=l j'=l 


(A2) 

(A3) 

(A4) 


and 


Ky/ 


iry,y — 
i-i'J-j' ~ 


iry,x 

^i-i', 3 - 3 ' 





y ~ dxdy 

[(x - a;,)2 + {y- y^Y^F 

(AS) 

(v - Vi)(v - Vi') 

[{x - XiY + {y- 

(A6) 

{x - Xi>){y - yj) 
l(x - X,Y + {y - 

(A7) 


The full expressions of force kernels used in this work can be found in Yen et ah (2012) and 
are summarized as follows for completeness: 
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where xi 1/2 ^i'+i /2 Vi yy —\/2 Vj and yu yji^ij2 yj- 

The corresponding expression of Ff- used in the particle-based method is the following: 
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where 


Vf - Vj 


or j ^ f 


1^0, otherwise. 

py;-^ ^ a.,Pll + 5lp^y + 5lp,%. (A16) 

Similarly, /Cg q = /Cg’o = 0 due to the odd symmetry with respect to the cell center. 


A.2. Polar coordinates 


Following Equation (24): 
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With the linear approximation in the surface density, Ffj with accuracy of 2nd-order can be 
approximated by: 
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Introducing some auxiliary symbols: 
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The full expressions of force kernels are: 
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B. Derivation of The Recursive Relations 


Without loss of generality, we set r = 1 to simplify the notation when deriving the 
recursive relation. The formula for r > a in Equation (53) can be recast as: 
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where Ti(q;) = a has been applied. Equation (52) then reads: 
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Integrate directly for A: = 0,1 in Equation (B3): 
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where we have used T 2 (a) = 2a'^ — 1 , T 4 (q;) = 8 a^ — 8 a^ + 1 and the relation 2 TmT„ = 
Tm+n + T\m-n\- Integration by part for the integral in Equation (B5) and apply the identity 
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Apply the change of variable a = cos(6*), where d = cos ^(a), and use the property 
Tn{cos9) = cos{n9) to the integral in Equation (B7); 
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where Un{ct) is the Chebyshev polynomial of the second kind of order n. From Equation (BIO) 
to Equation (Bll), we have used the relation: 
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The Chebyshev polynomial of the second kind can be expressed using the Chebyshev poly¬ 
nomial of the hrst kind through; 
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Therefore, Equation (Bll) can be expressed using the Chebyshev polynomial of the hrst 
kind: 
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From Equation (B17) to Equation (B18), the order of summation is exchanged. Combining 
Equations (B5)(B7)(B11)(B18), T can be expressed in terms of sin“^(Q;), \J\ — and the 
combination of Chebyshev polynomials: 
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Similar to the derivation for Equation (B19), it is straightforward to show that Equation (B4) 
has the following expression: 
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Substituting Equations (B19) and (B20) into Equation (B2), we obtain the recursive relations 
as shown by Equations (54) and (55). 


REFERENCES 

Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 








































- 28 - 

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University 
Press) 

Bracewell, R. N. 1999, The Fourier Transformation and its Applications, 3rd edn. (Me. Graw- 
Hill) 

Elmegreen, D. M., Elmegreen, B. G., Erroz-Ferrer, S., et al. 2014, ApJ, 780, 32 

Evans, L. G. 1991, Graduate Studies in Mathematics, Vol. 19, Partial Differential Equations 
(Providence, Rhode Island) 

Hockney, R. W., & Eastwood, J. W. 1988, Gomputer simulation using particles 

Inutsuka, S.-i., Machida, M. N., & Matsumoto, T. 2010, ApJ, 718, L58 

James, R. A. 1977, JGoPh, 25, 71 

Kim, W.-T., Seo, W.-Y., & Kim, Y. 2012, ApJ, 758, 14 

Lee, W.-K. 2014, ApJ, 792, 122 

Lee, W.-K., & Shu, F. H. 2012, ApJ, 756, 45 

Li, S., Buoni, M. J., & Li, H. 2009, ApJS, 181, 244 

Lin, L.-H., Wang, H.-H., Hsieh, P.-Y., et al. 2013, ApJ, 771, 8 

Mestel, L. 1963, MNRAS, 126, 553 

Miiller, T. W. A., Kley, W., & Mem, F. 2012, A&A, 541, A123 

Schulz, E. 2009, ApJ, 693, 1310 

Seo, W.-Y., & Kim, W.-T. 2014, ApJ, 792, 47 

Yen, G.-G. 2014, SJAM, 2014 

Yen, G.-G., Taam, R. E., Yeh, K. H.-G., & Jea, K. G. 2012, JGoPh, 231, 8246 
Zhang, H., Liu, H.-G., Zhou, J.-L., & Wittenmyer, R. A. 2014, RAA, 14, 433 
Zhang, H., Yuan, G., Lin, D. N. G., & Yen, D. G. G. 2008, ApJ, 676, 639 


This preprint was prepared with the AAS lAIpX macros v5.2. 



29 


Table 1: Order of convergence in terms of and L°° for the radial and azimuthal forces 

for those algorithms discussed in Figures 5 and 6. The numbers are extracted using the 
numerical errors obtained for Nd = 64,128, 256,512,1024. 


disk model 

error norm 

particle 

SIM 

particle-1- 

BM08 

^D2 


1.0 

1.5 

1.4 

1.0 


T OO 

1.0 

1.4 

1.1 

1.0 


Lr 

1.0 

1.9 

1.2 

1.0 


LI 

1.0 

2.0 

1.1 

1.0 


Ll 

1.0 

2.0 

1.2 

1.0 


Li 

V 

1.0 

2.1 

1.2 

1.0 


T OO 

l^r 

1.0 

1.9 

1.0 

1.0 


T OO 

1.0 

2.0 

1.1 

1.0 


Ll 

1.0 

1.9 

1.0 

1.0 


Ll 

1.0 

1.9 

1.1 

1.0 


Ll 

1.0 

1.9 

1.1 

1.0 


Ll 

1.0 

1.9 

1.1 

1.0 



Fig. 1.— Evaluation of /Cqo, /Cqo and /Cq’q at the cell center, as denoted by the 

asterisk symbols. The fan-shaped cell is covered with Cartesian cells. The forces are evalu¬ 
ated using the 2nd-order scheme described in Section 3.1 to avoid the singularity appearing 
in polar coordinates, (a) Surface density with a = 1 is used to evaluate /Cgg. (b) Surface 
density with unit slope in radial direction is used to evaluate rm/Cgo. (c) Surface density 
with unit slope in azimuthal direction is used to evaluate /Cq’q . The fan-shaped cell, which 
is characterized by Ax and A|/, should be spatially resolved using roughly 10 Cartesian cells 
to have a reasonable speed of convergence. 
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Fig. 2.— Comparisons of radial forces for the model with a = 0.25 and N = 128 in 
Cartesian coordinates. The solid line is the analytic solution, diamond symbol is the solution 
from the 2nd-order scheme, empty circle (particle+) is obtained using particle-based method 
with /Cq’o and /Cq’q correction and asterisk (particle) is the particle-based method without 
density slope correction, (a) Radial forces as a function of radius, (b) Relative error as a 
function of radius. 
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A^d 

Fig. 3.— The one norm error, L].^ as a function of cell number N for the (Td 2 model 
with a = 0.25 in Cartesian coordinates. The diamond symbol is the error from the 2nd- 
order scheme, empty circle (particle-|-) is obtained using particle-based method with density 
slope correction and asterisk (particle) is the particle-based method without density slope 
correction. The solid and the dashed lines indicate the slopes of Ist-order and 2nd-order 
convergence, respectively. The order of convergence htted for the last four data points are 
2.0, 1.1, 1.0 for the 2nd-order, the particle-|- and the particle methods, respectively. 
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Fig. 4.— The one norm and maximum norm errors as a function of cell number N for the 
aD 2 niodel with a = 0.006 in cylindrical coordinates. The diamonds and open circles are 
obtained from the particle-based method, while the asterisks and plus signs are from the 
original method proposed in Yen et ah (2012), i.e., = 1. The solid line indicates the 

slope of 2nd-order convergence. The order of convergence fitted for the last four data points 
are 1.8, 1.8, 2.0, 2.0 for the diamond, the open circle, the asterisk and the plus sign data, 
respectively. 







33 




Fig. 5.— The maximum norm errors as a function of cell number A'^d for the Udj rnodel with 
a = 0.25 in cylindrical coordinates. The center of 002 is placed at (xc = 0.5, r/c = 0.1). The 
maximum errors of the radial forces are shown in (a) and that of the azimuthal forces are 
shown in (b). The open circles and asterisks are obtained from the particle-based method 
without and with density slope correction, respectively. The diamonds are from the SIM 
proposed in Section 4 with iVtpz = 19. The inverse triangles are obtained using the method 
described in BM08. The red and blue lines indicate the slope of Ist-order and 2nd-order 
convergence, respectively. 
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Fig. 6.— The maximum norm errors as a function of cell number A'^d for the model with 
a = 0.25 in cylindrical coordinates. The center of is placed at (xc = 0.5, r/c = 0.1). The 
maximum errors of the radial forces are shown in (a) and that of the azimuthal forces are 
shown in (b). The open circles and asterisks are obtained from the particle-based method 
without and with density slope correction, respectively. The diamonds are from the SIM 
proposed in Section 4 with iVtpz = 19. The inverse triangles are obtained using the method 
described in BM08. The red and blue lines indicate the slope of Ist-order and 2nd-order 
convergence, respectively. 
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Fig. 7.— The maps of absolute error for the and disks discussed in Figures 5 and 
6. These maps are obtained using SIM proposed in Section 4. The left column is the results 
for the radial forces, while the right column is for the azimuthal forces. The top panel is for 
the aD 2 disk, while the bottom one is for the cr^g disk. 






















